A reliable benchmark of the last 640,000 years millennial climate variability

How often have past climates undergone abrupt transitions? While our understanding of millennial variability during the past 130,000 years is well established, with precise dates available, such information on previous climate cycles is limited. To address this question, we identified 196 abrupt transitions in the δ18O record of the well-dated Chinese composite speleothem for the last 640,000 years. These results correspond to abrupt changes in the strength of the East Asian Monsoon, which align with the Greenland stadials and interstadials observed in the North Atlantic region during the last 130,000 years before present. These precise dates of past abrupt climate changes constitute a reliable and necessary benchmark for Earth System models used to study future climate scenarios.

Abrupt transitions in the Earth System have become a major concern in recent decades, due to increasing evidence provided by paleoclimatic records 1 .These transitions can be classified into different categories, with Tipping Points (TPs) representing ultimate jumps from one climate state to another, with no possibility of returning to the initial conditions.In recent years, TPs have been studied with growing interest due to the impact of anthropogenic activities on the various components of the Earth system and their potential future evolution, leading to conditions that have not been observed over the last 2.6 million years.Recent analyses of current climate components, under different greenhouse gas emission scenarios, suggest that the Earth may already be approaching some of these tipping points, and that others are expected over the current century 2 .However, these projections are based solely on very recent data, dating back only a few hundred years.
The Earth system is considered to be reacting relatively abruptly to current anthropogenic forcing.However, the notion of abruptness remains ambiguous, as it refers to a timescale that is difficult to assess properly.Aware of this problem, the tipping elements listed by Lenton et al. 3 , which correspond to particular domains of the Earth system, are based on long-term observations under controlled conditions that make it possible to identify the associated tipping points.For example, there is evidence that if the rate of deforestation due to forest fires and climate change is not reduced, the Amazon rainforest will soon reach a tipping point towards a savanna state 4 .This would have a direct or indirect impact on regional and global climate systems, as well as on various other ecosystems 5 .Lenton et al. 3 have identified tipping elements that are closely linked to current climate change and therefore directly or indirectly related to anthropogenic forcing.However, interpretations of current abrupt climate changes must always refer to studies of abrupt climatic transitions evidenced in paleoclimate proxy records covering longer timescales with little or no human impact.The study of abrupt climate changes in the past has therefore become a new field of research in recent years 1 , exploring fluctuations that occurred in relatively short time intervals of several tens or, at most, hundreds of years, according to high-resolution records from the Greenland ice cores [6][7][8] .However, there is evidence that abrupt transitions can also be identified in the deeper time from lower resolution records, which still reveal changes or transitions that have had a considerable impact on the global dynamics of the Earth system [9][10][11] .

Abrupt transitions in past climates
Recent studies of past climate changes have underlined the importance of studying abrupt transitions in the past, while recognizing the challenges associated with data quality, including timescale accuracy and quantification of associated uncertainties 12 .Studying past abrupt transitions and the mechanisms involved requires the best

The abrupt climate transitions of the past 640 kyrs
To assess abrupt climate transitions of the past, well-dated, high-resolution records and robust statistical approaches are needed.Greenland ice cores provide such records, capturing abrupt transitions between cold and warm conditions known as Greenland Interstadials (GI) or Dangaard-Oeschger (DO) events, with a 20-year resolution that enables precise dating of these transitions 7 .Dates have also been proposed for transitions to colder conditions, known as Greenland Stadials (GS).However, the timescale of these abrupt transitions is limited to the last climate cycle, which hinders the evaluation of Earth System models over a longer time interval.To overcome this limitation, Greenland δ 18 O variations were reconstructed from EPICA δ 18 O measurements using the thermal bipolar seesaw model of Stocker and Johnsen 32 .This method enabled the detection of abrupt, GI-like warmings over the past 800,000 years with a time resolution of 50 years 33 .However, the reconstructed abrupt warming events do not follow the decadal timescale identified in Greenland for the last climate cycle.In-depth testing of Earth System models may therefore require even finer resolution records.Another high-resolution record that can be used for this purpose is the δ 18 O record of the Chinese composite speleothem (CS), which covers seven climate cycles over the past 640,000 years 34 .Its high resolution and precise dating make it one of the best paleoclimate records for detecting abrupt transitions over several climate cycles and providing a reliable benchmark for the study of past abrupt climate changes.
As paleoclimate records differ in origin, duration, and periodicity, it is essential to have an objective, automated methodology for identifying and comparing abrupt transitions.Among several methods used to detect such transitions, a new approach based on the nonparametric Kolmogorov-Smirnov (KS) test (e.g. 35) has recently been developed and has shown promising results 36 .The KS test is applied to compare two samples drawn from a time series, one before and one after a potential jump.The KS statistic measures the difference between the empirical distribution functions of the two samples, enabling discontinuities in the time series to be identified.To refine the results and pinpoint significant transitions, the KS test is augmented by additional criteria, including a variable window size and a minimum rate-of-change threshold.Window size is a critical factor affecting the identification of transitions (see 14 for more details).In addition to the KS test, long-term trends in maxima and minima can also be used to establish key transitions, such as Stadial-Interstadial boundaries in the Greenland ice cores.KS method's parameters are optimized using receiver operating characteristic (ROC) analysis.Using this accurate and robust approach, it can be noted that transitions in paleoclimate records are often not correctly identified in the literature, possibly due to variable data quality and dating methods.Although the KS test is sharper and can find precise transition dates, recurrence quantification analysis (RQA) is better suited to identifying particularly important transitions that correspond to changes in dynamics.To quantitatively assess these major transitions, recognized as distinctive patterns in Recurrence Plots (RP), we employ an RQA measure known as Recurrence Rate (RR).The minima of RR, identified through prominence analysis 36 and marked by pink crosses in Fig. 3, signify substantial shifts in the system and are thus of particular interest to us.
The use of a variable window size in the KS test can enable more accurate identification of transitions occurring over different time scales in paleoclimate records, revealing transitions that may have been missed in other analyses.In the case of the CS δ 18 O record, the KS test was applied using two different window size ranges (0.4-4 kyr and 0.6-4 kyr) to optimize transition identification (Fig. 1, Table 1).This approach detected a total of 144 similar transitions over the last 640 kyrs, as well as 27 transitions having different dates between the two analyses.While reducing the window size diminishes the statistical significance, it enables a more precise examination of the time series.As a result, the window size of 0.4-4 kyrs allowed the identification of additional 25 transitions to both moist (11) and drier (14) conditions.This improvement is not just cosmetic, as it has also enabled the detection of events found in other paleoclimate records, such as sub-events in the Greenland NGRIP δ 18 3).
for the last climate cycle 7,10,36 .Specifically, the KS test of the CS δ 18 O record with the 0.4-4 kyrs window detected 58 NGRIP events and sub-events, including the GSs and GIs (Table 2).To compare detected CS transitions to strong monsoon conditions with NGRIP GIs, we followed the labeling of Chinese interstadials detected in the last 130,000 years as "A" by Cheng et al. 19 and Wang et al. 21.A close correspondence appears between A24 and GI-24.2 at the base and A1 and GI-1e at the top of the records (see Fig. 2, Table 2).This correspondence is also observed for older GIs found in the reconstructed Greenland δ 18 O 33 .Following the original numbering of the GIs 37 , labelled as DO events, the KS test detected 23 out of 24 of such events, with only GI-9 or A9 missing.The difficulty in identifying transitions between GI-24.2 and the base of the last climate cycle, between 135,550 years and 128,550 years BP, may be due to the complexity of the correspondence between NGRIP and CS records, as previously shown by Rousseau et al. 10 .This difficulty could also be due to the fact that these paleoclimate records from different regions may respond differently to global climate changes, leading to regional variability in the timing and magnitude of abrupt transitions.
In the penultimate cycle (243-130 ka BP), the Chinese interstadials were labeled "B" by Wang et al. 21.However, the KS test did not detect the most recent moist events B6 to B1 (144-131 ka BP) while identifying most of the oldest (Table 1, Extended Table 1).Comparing abrupt transitions to moister events with reconstructed GIs from the Barker et al. 33 study, there is very little correspondence between the two records, with only 48 events, out of the 103 identified by Barker et al. 33 , having relatively close dates (see Extended Table 1).As previously mentioned, the reconstruction of Barker et al. 33 does not reproduce the decadal time scale of the NGRIP record 7 for the last climate cycle, leading to such result.When identifying strong monsoon events using the climate cycle boundaries defined by Lisiecki and Raymo 38 , the last two cycles have the highest number of abrupt transitions, with 21 and 28 identified for the penultimate and last cycle, respectively.In contrast, the four previous cycles, from 621 to 243 ka BP, never exceeded 11 abrupt transitions each, with 11, 9, 10, and 10 events respectively.This difference could be linked to the duration of the cycles, as only the last two cycles exceeded 110 kyrs.In addition, temporal resolution varies with climate cycles and could be considered to have a potential impact on the detection of abrupt transitions in older climate cycles.Our compilation of the composite speleothem record, using the Lisiecki and Raymo boundaries 38 shows that this does not appear to be the case (see Extended Data Table 1).
Recurrence quantification analysis carried out on the CS δ 18 O record exhibits a drift topology characteristic of a monotonic trend in variability, which corresponds to a clear dominant monsoon signature of 20 kyrs (Fig. 3).The analysis of the recurrence rates identifies 34 significant minima with a recurrence rate prominence higher than 0.5, already which correspond to transitions also identified by the KS test (Extended Data Table 2).Two main groups can be observed, separated by a main transition at 333 kyrs, with the younger transitions being more abrupt transitions than the older ones.Interestingly, the last 336 kyrs represent the time interval during which the greatest amplitude between interglacial and glacial is observed in the global δ 18 O record 38 .The analysis of the recurrence rates detects a key transition at 226.5 ka BP, corresponding to Chinese event B24, which appears to be major in terms of variability dynamics (see Fig. 3, Extended Table 2).The mechanism behind such dichotomy over the past 640 kyrs needs to be studied in greater depth using Earth system models.The KS test is also a useful tool for determining the start and end dates of various Terminations (T) that mark abrupt transitions between glacial and interglacial stages.The precise dates of these Terminations are shown in Table 3, indicating an average duration of 6.2 kyrs ± 1.9 kyrs.However, Cheng et al. 34 noted that these events are complex and often include one or more events corresponding to intervals of strong monsoon, surrounded by periods of weak monsoon known as "Weak Monsoon intervals" (WMI).Our analysis indicates that in the CS record, all seven identified Terminations begin with an abrupt drop in the δ 18 O signal, corresponding to the occurrence of a WMI.While five Terminations (TI, TII, TIII, TIV, TVII) show abrupt summer monsoon strengthening (high δ 18 O), characterized by either a very short or longer interstadial-like interval, two Terminations (TV and TVI) do not show such a complex structure (see Fig. 1, Table 3).This difference can be explained by the fact that these two Terminations are the shortest and correspond to low precession maxima and low insolation at 65° N (Berger 39 ), compared to the other Terminations (see 34 , Fig. 1).Earth System models should nevertheless study in detail such key issue.
The "2 kyr shift", discussed by Cheng et al. 34 and corresponding to a strong summer monsoon in East Asia during the late Holocene, was not detected by the KS-test.However, the test did identify a weak monsoon transition at 3.5 ka BP.This corresponds to a global cooling event 40,41 , which led to a significant advance of glaciers in Central Asia, the Southern Hemisphere, Northern America, Scandinavia, and the Alps [42][43][44][45] .It is also associated with tropical aridity in East Africa, South America and the Caribbean 46 , as well as with major atmospheric changes, such as the strengthening of westerlies in the North Atlantic Ocean 47 and the increasing strength of the Siberian High 48 , which constrains the East Asian Winter Monsoon 49 .

Implications for the analysis of past climates
What factors can be attributed to the abrupt changes detected in the composite Chinese Speleothem (CS) δ 18 O record?While the millennial variability reconstructed by Barker et al. 33 provided possible dates for abrupt warmings during the past 800,000 years, our analysis provides precise dates of abrupt transitions from the CS δ 18 O record using a robust statistical method.These transitions are linked to East Asian summer monsoon variability, with higher δ 18 O values corresponding to weak GS-equivalent intervals, and lower δ 18 O values corresponding to strong DO events-like equivalent intervals.Furthermore, the close correspondence between abrupt transitions detected in both the CS δ 18 O and Greenland NGRIP records over the last climate cycle (last 130 kyrs) leads us to propose that all precise abrupt transitions detected by the CS serve as a benchmark for subsequent analyses of the 640 kyrs millennial variability by Earth System models to better predict future climate change scenarios.
Using both the augmented Kolmogorov-Smirnov (KS) test and Recurrence Quantification Analysis (RQA), our study aims to provide a comprehensive analysis of abrupt climate transitions in the composite Chinese Speleothem record.The augmented KS test is used to identify significant changes in the data, while RQA and the associated Recurrence Plot (RP) provide further insights into the recurrence patterns and system dynamics.
The first approach to identifying abrupt transitions in our datasets is to use the augmented KS test as described in Bagniewski et al. 36 18 O record with the KS test using the 0.4-4 kyr window range; similar events published for the NGRIP δ 18 O record by Rasmussen et al. 7 with their corresponding labels; labels used by Cheng et al. 34 for the published strong monsoon intervals; abrupt warmings reconstructed by Barker et al. 33 using the EPICA δ 18 O record and the corresponding labels; comparison between the CS abrupt transitions detected with the KS test and the reconstructed Greenland interstadials from Barker et al. 33 .GS Greenland Stadial, GI Greenland Interstadial, CI Chinese label.Same conventions than in Table 1.
time t, (X t ), using varying window lengths.The KS statistic is calculated for all window lengths to detect abrupt transitions, with values above 0.7 considered significant.Transition detection is subsequently refined using a minimum rate-of-change threshold.The analysis begins by identifying transitions using the longest window, which has the highest sample size and therefore the greatest statistical significance.The method then integrates transitions detected using shorter windows to capture transitions on shorter time scales.To identify transitions between dominant climate modes, such as the boundaries between weak and strong monsoons in the Chinese speleothem record, we use a running window to determine the upper and lower values of the time series and identify transitions that correspond to a shift from one mode to another.Two window length ranges, 0.4-4 kyrs   www.nature.com/scientificreports/RPs are a powerful tool for identifying recurring patterns in time series data, including paleoclimate records.The RP for a time series takes the form of a square matrix whose two axes represent time.A dot is inscribed at a position (i, j) in the matrix when |x i − x j |< ε, ε being the recurrence threshold.The resulting matrix provides a visual representation of recurring patterns in the time series.In addition to the visual interpretation of RP, the RQA is used to objectively quantify the recurrence patterns.Eckmann et al. 50distinguished between large-scale typology and small-scale texture when interpreting a square matrix of dots created by recurrence plots (RP).The most intriguing RP typology concerns recurring patterns that lack periodicity and are challenging to detect through conventional spectral analysis techniques.Marwan et al. 51 reviewed the various methods to objectively quantify these visual RP typologies, collectively known as recurrence quantification analysis (RQA).One of the parameters derived from RQA is the Recurrence Rate (RR), which describes the probability of recurrence of system states in a given time interval.In this study, the RR is calculated using a sliding window of 4 kyrs and plotted alongside the RP.The minima of the RR (Extended Table 2), identified as abrupt transitions, are selected according to their respective prominence following the approach by Bagniewski et al. 36 .
The combination of the two methods provides a robust and comprehensive analytical capability.It has been successfully demonstrated and applied to a variety of paleoclimate records spanning the last 130,000 years to the last 66 million years from different domains (glacial, marine and continental) [9][10][11]14,36,52 .
Table 3. Dates of the climate Terminations detected in the CS δ 18 O record by the KS-test.From left to right: termination number according to the marine isotope stratigraphy; start and end dates of the Terminations; duration of the Terminations, including their mean duration, standard deviation, and maximum and minimum values.Four additional intervals indicated by Cheng et al. 34

Figure 1 .
Figure 1.Detection of the abrupt transitions towards weak and strong East Asian Monsoon intervals.δ 18 O variations during (A) the 0-220 kyrs interval; (B) the 220-440 kyrs interval; and (C) the 440-640 kyrs interval with the upper and lower panels showing transitions detected with the KS test applied with the 0.6-4 kyr and 0.4-4 kyr window size range, respectively.Blue lines mark the transitions towards weak monsoon, red lines mark the transitions towards strong monsoon.Indication of the different Terminations, T1 to TVII with mention of TIII-A (#) and TVII-A ( §) (see Table3).

Figure 3 .
Figure 3. Recurrence Quantification Analysis (RQA) of the Chinese Speleothem δ 18 O composite record.(A) δ 18 O record over the last 640 kyrs BP, (B) recurrence plot (RP) of the δ 18 O time series, (C) recurrence rate (RR) versus time with minima marked by pink crosses selected according to their respective prominence (see Extended Table2).
O record

Table 2 .
. The two-sample KS test involves comparing values on either side of a value of a proxy X at Comparison of the abrupt transitions detected for the last climate cycle (All dates in ky).From left to right: transitions detected for the Chinese Composite Speleothem δ Comparison of the detected abrupt climate transitions in both NGRIP and the Chinese Speleothem δ 18 O records over the last climate cycle (see Table 2).A, during the 46-8 kyrs interval; B, the 84-46 kyrs interval; and C, the 122-84 kyrs interval.Upper panel NGRIP and lower panel Chinese Speleothem.